function [x,normk,it,nrm] = redblack2D(x,maxiter)
% function [x,normk,it,nrm] = redblack2D(x,maxiter)
%    x     = [jout' keff]'
%    normk = difference in consecutive keff's
%    it    = number of outer iterations
%    nrm   = norm of nonlinear residual, f([x 1]), where the nonlinear
%            function is the one used for the Newton solver
% j. roberts 2/22/2010

% global variables set:
%   order = discrete legendre order for the response function expansion
%   numg  = number of energy groups
%   numel = number of elements
global order numg numel

% 1 element case
if numel == 1
mm = [ 1 0 0 0 
       0 0 0 0
       0 0 0 0 
       0 0 0 0 ];                                % mm maps faces to faces
else
% 2 element case
mm = [ 1 0 0 0 0 0 0 0
       0 0 0 0 2 0 0 0
       0 0 1 0 0 0 0 0
       0 0 0 1 0 0 0 0
       0 2 0 0 0 0 0 0
       0 0 0 0 0 1 0 0
       0 0 0 0 0 0 1 0
       0 0 0 0 0 0 0 1];      
end

% tmp goes where reflection occurs
tmp  = diag( ones(1,order+1).*(-1).^(1+(1:order+1)) );
tmp  = kron( tmp, eye(numg) );
% tmp goes where transmission occurs
tmp2 = diag( ones(1,order+1).*(1).^(1+(1:order+1)) );
tmp2 = kron( tmp2, eye(numg) );
% M combines everything
M   = kron(mm==1, tmp )+kron(mm==2, tmp2);
M   = sparse(M);
N   = length(M(1,:));  % should be 8 for 0th, 16 for 1st, etc.
j1  = M*x(1:end-1, 1);
k1  = x(end, 1);
jred = zeros(N,1); jblk = jred;
jred(1:2:end)  = j1(1:2:end);
jblk(2:2:end) = j1(2:2:end);
k0 = 0; 
it = 0;
jwarn=0;
jjit=0;
maxjit = 6000;
jold=zeros(N,1); nrmm = norm( feval('respfct2Dexample',[j1' k1 1.0]') );
while nrmm > 1e-10 && it < maxiter%norm(k1-k0) > 1e-9 && it < 2000
    disp(['order/it: ',num2str(order),' ',num2str(it)])
    % get response matrices for current k value
    j0 = zeros(N,1);
    [R1 F1 A1 L1] = resp2dexample( k1, order, 1  );
    R1 = sparse(R1);
    if (numel==1)
        R = R1; F = F1; A = A1; L = L1;
    else
        % if two element
        [R2 F2 A2 L2] = resp2dexample( k1, order, 2  );
        R2 = sparse(R2);
        R = sparse( [R1 zeros(size(R1)); zeros(size(R2)) R2] );
        F = [F1 F2];
        A = [A1 A2];
        L = [L1 L2];
    end

    % else
    % R=R1;F=F1;A=A1;
    jit = 0;
    jold = j1;
    while norm(j1-j0) > 1e-12 && jit < maxjit
        sumj = sqrt(j1'*j1);
        j0      = j1;
        j1      = M*R*j1/sumj;
%         jblk    = M*R*jred/sumj;
%         jred    = M*R*jblk/sumj;
%         j0      = j1;
%         j1      = jblk + jred;
        jit = jit + 1;
    end
    if (jit==maxjit)
        jjit=jjit+1;
        jwarn=max(1,jwarn);
    end
    k0 = k1;
    % keep in mind j's are INCOMING currents
    jout = R*j1; jout = jout/sqrt(jout'*jout);
    gain = F*j1;  
    LL=0; LR = 0; LB = 0; LT = 0;
    if numel == 1
    i1=1; i2 = N/4; i3 = N/2; i4 = 3*N/4; i5=N;
    LL = L(1,:)*j1(i1:i2)+L(2,:)*j1(i2+1:i3)+...  % left side leakage
         L(3,:)*j1(i3+1:i4)+L(4,:)*j1(i4+1:i5);
    LB = L(1,:)*j1(i3+1:i4)+L(2,:)*j1(i4+1:i5)+...
         L(4,:)*j1(i1:i2)+L(3,:)*j1(i2+1:i3);
    LR = L(2,:)*j1(i1:i2)+L(1,:)*j1(i2+1:i3)+...  % right side leakage
         L(4,:)*j1(i3+1:i4)+L(3,:)*j1(i4+1:i5);
    LT = L(2,:)*j1(i3+1:i4)+L(1,:)*j1(i4+1:i5)+...
         L(3,:)*j1(i1:i2)+L(4,:)*j1(i2+1:i3);
    else
        % NN=N/2;
        % i1=1; i2 = NN/4; i3 = NN/2; i4 = 3*NN/4; i5=NN;
        % LL = L1(1,:)*j1(i1:i2)+L1(2,:)*j1(i2+1:i3)+...  % left side leakage
        %      L1(3,:)*j1(i3+1:i4)+L2(4,:)*j1(i4+1:i5);
    end
    loss = A*j1 + abs(1-mm(1,1))*LL + abs(1-mm(2,2))*LR +...
                  abs(1-mm(3,3))*LB + abs(1-mm(4,4))*LT;
    k1 = gain/loss
    if (k1 < 0.001)
        disp('problem')
    end
    it = it + 1;
    nrm(it) = norm( feval('respfct2Dexample',[jout' k1 1.0]') );
    nrmm = nrm(it);
    normk(it) = norm(k1-k0);
end
nrm = [1 nrm];
if(jwarn>0)
disp('warning: at least one inner did not converge')
disp(['# of kits: ',num2str(it),' # of jjits ',num2str(jjit)])
end
x = [jout' k1]';

